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Lagrangian technique associated with two different time schemes: explicit and sem!- 
implicit. Both schemes produce realistic fronts after approximately 40 hours of model 
integration. The semi-Lagrangian semi-implicit scheme ıs more successful in handling 
the sharp gradients associated with the front. Also, the semi-Lagrangian semi-implicit 
equations are integrated with time steps as long as 3600 sec. producing solutions with 
relatively smallerrors. This indicates that this numerical scheme is appropriate for use 


in mesoscale regional models. 
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|. INTRODUCTION 


The development of fast computers has been of great importance for the 
atmospheric sciences. With them, scientists have been able to obtain a better 
understanding of the physical processes tn the atmosphere and operational weather 
services have produced more accurate forecasts. But as pointed out by Robert (1981), 
scientists can also give their contribution by developing numerical algorithms that are 
computationally less expensive and at the same time accurate. This interaction between 
science and technology will allow the computational effort to be used in numerical 
models that will have higher resolution and that will be more accurate, as well, since 
more realistic and complex physical processes would be represented in these models. 

One of the atmospheric processes that has been under active study is the formation 
of fronts. In a study of non-geostrophic frontogenesis Williams(1972) reproduced 
numerically a realistic front, demonstrating the importance of nonlinear effects in this 
dynamical process. Since frontogenesis is a phenomenon that ts characterized by the 
development of sharp gradients of the physical parameters it would be interesting to 
develop a numerical model based on techniques that can handle adequately the scale 
collapse problem. Kuo and Williams (1990) studied the use of different numerical 
techniques in the solution of a simple scale collapse problem. They concluded that the 
semi-Lagrangian method, developed by Sawyer (1963), would be adequate for use in 
numerical models where strong gradients were present. 

The semi-Lagrangian technique has been used mostly in synoptic scale numerical 


models. Robert (1981) used this technique with the semi-implicit scheme (Robert et al.. 


1972) 1n a primitive equation model and obtained perfectly stable solutions using time 
steps 25 times larger than those that would be allowed by the Courant-Friedrichs-Lewy 
(CFL) stability criterion. Pudykiewicz and Staniforth (1984) examined the use of the 
scheme in several different conditions and concluded that besides its good stability 
properties, it is also accurate and flexible. 

The objective of the present study is to perform experiments using the semi- 
Lagrangian technique in a frontogenesis model and verify its suitability for use in the 
representation of atmospheric processes where sharp gradients are present. 

The semi-Lagrangian technique is introduced tn the next chapter. The basic model 
equations are developed in chapter III. The numerical procedure employed for solving 
the equations ıs presented in chapter IV. In chapter V the experiments and results are 


discussed, followed by the summary and conclusions in chapter VI. 


No 


Il. THE SEMI-LAGRANGIAN SCHEME 


This model uses a semi-L agrangian algorithm similar to the one described by 
Robert (1981). 

The three time level semi-Lagrangian technique consists of evaluating the total 
derivative of a general dependent variable Q(F(t),t) following the trajectory of the 
fluid particle, using the approximation 


dQ _ Q(r(t+Adt),t+At)-Q(r(t- Ad),t- At) (2.1) 
dt 2At 


where F(t) represents the position vector of the fluid particle at time t£. For sake of 
simplicity, constrain the problem to two dimensions, y and Z. In this case the position 


vector will be 


F(t) = y(Dj+2z(DŘ. ae 


The location of the particle at the forecast time (t + 4 t)is chosen to be coincident 


with a grid point. Let the displacements of the particle along the y and Z directions 


during one time step be represented by a and b. respectively. Ifthe grid point position 


is represented by 


r(t+ At) = OZ) (2.3) 


the approximation (2.1!) can be written as 


dQ QO p% tt Ad - Q(y, - 24,2, -2b,t- At) (2.4) 
dt 2At | 


The displacements a and b can be calculated iteratively using 


a™! -~ At.v(y,-a*,7,-b",t) l (2.58) 


+] 2 
pro = At.w(y,-a",z,-6",t), (2.5b) 


where the superscript o represents the n” iteration and v and w represent the 
horizontal and vertical velocities, respectively. The initial estimates of aand bin the 


iterative process are defined by 


0 = I Se 
a = At.V(V;Zp 8), (2 Sc) 


b'— At.w(3,,2,:8) ee 


[n Pudykiewicz and Staniforth (1984) 1t can be found that the necessary condition for 


convergence of the iterative procedure described in (2.5) is 
BB) e 
ov) log |oy| jo 


Furthermore. Kuo and Williams (1990) point out that no more than three Iterations are 


Ar. max 





necessary for convergence to be attained. Robert (1981) compared the use of two and 
four tterations in his primitive equation model and the results obtained were quite 
similar. 

In order to determine the values of the variable Q at time t-At it is necessary to 
usean interpolation scheme. Ín this model, bicubic spline interpolations (de Boor. 1962) 
were used with the explicit formulations of the semi-Lagrangian scheme, using the 
algorithm described by Marchuck (1982). In the semi-ımplicit formulation only one- 
dimensional cubic splines in the y direction are necessary to solve the problem: as will 


be shown in IV.C.2. 


KA 


HI. BASIC EQUATIONS 


This model uses the hydrostatic Boussinesq equations. which assume 
incompressibility of the atmosphere. Friction and heating are neglected. It is also 
assumed that the atmosphere is bounded above by a rigid lid. Periodic boundary 
conditions are used in the horizontal. 


The Boussinesq equations can be written as 


du dq i 
L ee : (3.1) 
a Ca? 

O ee 3.2) 


= 0, (3.3) 


- 0, (3.4) 





-_—, (3) 
Oz sj 
where 
E 3.6) 
C, 
0 -T(Ê) -6,, (3.7) 
Po 
E PY 3 
p= COM) +87. (3.8) 
Po 


This model assumes that the Coriolis parameter fis constant, and that the bottom 
topography is flat. The variations in the x direction, along the front, are neglected, 
except for the q field and the basic deformation field. 


With the assumptions above, equation (3.3) reduces to 


ov | Ow (3.9) 


Ea 


dy 


for departures from the deformation wind field. 


With the use of a rigid lid assumption the pressure ts not known at any level. and 
this will require some manipulation of the equations. 
Integrate the hydrostatic equation (3.5) from z = Oto an arbitrary level z and 


obtain 


(3.10) 


$ 
| 
o w N 
DP 
R- 
+ 
A 


The constant of integration, C, can be determined by taking the vertical average of 


equation (3.10) and subtracting from equation (3.10) which gives 


Z Z 
ọ = É [ode-<É [0de>+<p>, (3.11) 
907 00. 


where < > denotes the vertical average operator 


H 
<F> = =$ Faz. (3.12) 
Hy 


The constant of integration has been eliminated, but now an expression for <@> must 


be found. 


Apply the vertical average operator to the continuity equation (3.9) and use the 


vertical boundary conditions 


w = 0 at 2=0 (3. 13a) 
and 
w-0 a z-H a) 
which vields 
E aj (3.14) 





Next. expand the total derivative in equation (3.2), multiply (3.9) by v and add 


to (3.2) to obtain 


- -— -fu. (313) 





Take the vertical average of (3.15), using the vertical boundary conditions (3.13), 


which gives 


JEW e o<yy> E d<@> 


ot dy dy 


~f<u>. (3.16) 





Differentiating (3.16) with respect to vand using equation (3.14) yields 


<g> _ View - (3.17) 


dy ey" oy 


With appropriate boundary conditions, a solution for <p> can be found. 

In order to solve equation (3.1), an expression for = must be obtained. 
Differentiate (3.11) with respect to x and use the assumption that ĝis not a function of 
x which gives 

ea. jew 
Ox a 


(3.18) 





An expression for e must be obtained. Take the domain average of (3.2) 


(<2 >) z k- Po} afu), (3.19) 


where f( )! denotes the horizontal average operator 


Fdv. (3.20) 


—— nin 


(Fy - 1 
L 


wie 


Expand the first term of (3.19) and use the vertical boundary conditions (3.13) which 


yields 


dv, H< v> _ O<v> 


<>; 


dt ot ot 


(3.21) 





Using periodic boundary conditions in the y direction, it can be shown that the 
second term in (3.19) is zero. With the assumption of constant 4 and using (3.21), 


equation (3.19) becomes 





Cate - -fi<u>}. (3.22) 
ot 


Likewise. take the domain average of (3.1) and use the vertical and horizontal boundary 


conditions which yields 


dikur) MA A< p >) a (3.23) 
ot Ox 


Equations (3.22) and (3.23) show that if {<u>} is set initially equal to zero, <v> 


. E dé <p> , 
will remain constant in time and so will dear Furthermore, if <v> is also set to 


zero initially, ser will be identically zero for all time. 





. aa Bid | 
In this study it will be further assumed that =e Is constant in the cross frontal 


direction (y), which allows to write 


ANEP _ <p> (324) 


Ox Ox 


With the assumptions above the following relation is obtained 


o . aa a. (3.25) 
Ax Ax 


for initial conditions 


(3.26) 


| 
© 


{<u>} = <y> 


The dependent variables will be expressed as a combination of a basic and a 


disturbance part as follows: 


u(x,y,Z,t) - U,(x,y)+u'(y,z,0, (3.27) 
v(x,y,2,0) - V,(x,y)+v'/(y,2,0) (3.28) 
(3.29) 


W(V,Z, 8) oz w’(y,Z,0), 


Bie.) = Cy) (3 30) 


P (x,y,z) - (x,y) +p (Y,D, (331) 


where the quantities with primes represent the disturbance part of the variables. 
The basic wind field that drives the frontogenesis is the nondivergent and 


irrotational deformation field given by 


D 
U,(x,y) - - —* sinh(px) sin(py), (3.32) 
u 


Da 
V (2,0) Pon aaa cosh ( p x) cos(H y). (3.33) 
u 


This wind deformation field was originally used by Stone (1966) in his study of quasi- 
geostrophic frontogenesis. 

It has been assumed that the perturbations are independent of x If the 
expressions for u and v (3.27 and 3.28 respectively) are substituted into the basic 
equations, o and = will be functions of x. In this study, where the main interest 1s 
in examining cross front variations, the equations will be evaluated at x = 0. This will 


make the mathematical formulations compatible with the previous assumptions and it 


is expected that the results will be physically consistent. 


Therefore, the basic wind deformation field will become 


U (0,3) - 0, (3.34) 
V (0, y) = 2% ita (3.35) 
u 


The basic geopotential field @ is chosen such that it will be in geostrophic balance 


with the wind deformation field as follows: 


fU, - - ge (3.36) 


oy 


Vv, == ID (3.37) 


Substitution of (3.36) and (3.37) into (3.27) and (3.28) respectively, and use of the 
expressions obtained for the dependent variables in equations (3.1)-(3.4) yields the 


following set: 


du 
AT (3.38) 
dt J 
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OV 
a LO py (vV) 4, (3.39) 
dt oy dy 
wv, ow 0, (3.40) 
dy & 
2% q (3.41) 
dt 


where the primes were dropped from the disturbance variables. 


Finally, an easier way to obtain = will be described. The dependent variable g 


can be expressed as 


P = Pot, ie) 
where 
rgo 
P, = jse (3.43) 
o Vo 


Differentiate (3.42) with respect to v. to obtain 


op - ka g 8 (3.44) 
Oy Y dy D 


[he value of g,1s not known, by virtue of the rigid lid assumption. In practice v 
will be predicted using equation (3.44) with Evo assumed initially equal to zero. After 
that, the v field will be adjusted in order to satisfy < v > = 0. This procedure will be 
equivalent to solving equations (3.11) and (3.17) for q, differentiating the values 
obtained with respect to y, and substituting into (3.39). Equations (3.38)-(3.41) and 
(3.44) give a complete set that can be integrated in time in order to predict the values of 


the dependent variables. 


IV. NUMERICAL SOLUTION 


A. INITIAL CONDITIONS 


The initial conditions used in this model are similar to the ones used by Williams 


et al. (1991) in their study of effects of topography on fronts. 


The initial temperature field is given by 





YX, SE [z 71 -E 


where 


e - 0, ZEZ 


, Dale 


(4.1) 


(4.2) 


(4.3) 


In the expression above ——* is a constant, A is the amplitude of the temperature 


disturbance and z, is the height of the upper limit where the temperature disturbance is 


present. This initial temperature field will confine the frontogenesis to the lower layers 


of the atmosphere, which will make the results more realistic. 


The initial y component field is chosen to be in thermal wind balance with the 


temperature field and is given by 


SpA 


2 fO, sin(py) . A 


u(y, z,0) pa 


where 








À ~ 0, Zaz 


and the condition 


u(y,Z,,9) = 0 


has been used. 
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> 252, 


(4.4) 


(4.5) 


(4.6) 


(4.7) 


The initial v and w fields are obtained from the solution of the quasi-geostrophic 


circulation equation (Williams, 1972) 


§ om y, ey -28 Ya 98 (4.8) 
aor Oz dy? dz? o gy gy 
where 
ON (4.9) 
Oz 
w owed (4.10) 


The vertical boundary conditions for the solution of (4.8) are set to 


w(z - 0) - w(z- H) - 0 Ben 


It can be shown, from the domain averaging of equation (4.2) that {< u >} will 
be zero initially. Also, the condition (4.11) will allow (3.25) and (3.26) to be satisfied in 


order for the numerical solution to be consistent with the assumptions of the model. 


B. THE SPATIAL GRID 

The numerical calculations are performed on a staggered grid, both in the 
horizontal and in the vertical directions. In the horizontal the variables are staggered 
following the scheme B described in Arakawa and Lamb(1977). This scheme, as shown 
in Haltiner and Williams (1980), allows a better representation of geostrophic 
adjustment processes. The variables are also staggered in the vertical. Pielke (1984) 
points out that the staggering of the u and w components in the vertical gives better 
solutions for the vertical velocity from the continuity equation. 

The horizontal grid uses Í grid points separated by a constant grid space A y. The 
vertical grid has K grid points and also uses a constant grid spacing, referred to as Az. 


Figure | shows the arrangement of the variables on the grid. 


C. FORMULATION OF THE MODELS 


l. Semi-Lagrangian Explicit Scheme 
In the semi-Lagrangian explicit formulation of the present model, the total 
derivatives will be approximated using the technique described in chapter II. The 
remaining terms of the equations are evaluated at the positions of the fluid particles at 
time level ¢ Introducing those approximations into the prognostic equations (3.38), 


(3.39) and (3.41) the following set is obtained, after dropping subscripts /and y: 


U(V,Z, t-At)-u(v-2a,2-2b,t-At) 


- f.v(y-a,z-b,t), (4.12) 
Al fo ) 


O w z(k+1) 


Z(k-3) 


y (j-3) y(j+3) y (j+1) 





Figure 1: Staggered spatial grid 
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v(y,2,1+ At) - v(y-2a,2-2b,t-At) _ 
2At 


E A D-az-b0 -f.u(y-a,z-b,t) 


V 
~ [v(y-a,z- 5,0) a V ¿O -a)) 0-0, (4.13) 


OC y,z,t+ At)-6(y-2a,z-2b, t- At) _ 


0. (4.14) 
2At 


The values of wand q are obtained from the diagnostic equations (3.40) and (3.45), 


respectively. 


2. Semi-Lagrangian Semi-Implicit Scheme 
The semi-Lagrangian semi-implicit scheme was introduced in this model 
following the development described by Robert et al.(1985). The scheme consists of 
separating the temperature into a basic state, dependent only on the vertical coordinate 
Z, and a disturbance part. Next, an implicit treatment is given to those terms related to 
gravity waves, namely: horizontal pressure gradient, divergence and vertical motion. 
To obtain the system of equations in the semi-implicit formulation take the 


advective terms out of the total derivatives in equations (3.38), (3.39) and (3.44); 
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d.u 
= du (4 15) 





OV 
——— +W = - + -fu-(v+V,)—*, (4.16) 


al ti w 00 (4.17) 





where the terms with subscript H represent horizontal components of the total 
derivative. 
Replace the temperature and geopotential fields by basic state and 


disturbance parts as follows: 


0(v,2,0 - 0'(2)40'(v,z,8), (4.18) 


(9,20 — 9'(2)40'(D,Z,0, (4.19) 


where the superscripts (*) represent the basic state parts of the variables. 
Now introduce (4.18) and (4.19) into equations (3.5), (3.40) and (4.15)44.17), 


and take the time average of those terms related to gravity waves, which yields 


ty 
ad 


dv aa 
tii wi a (vv), 
at OZ oy 

2) a 

dy a 














+ w— + = Q, 
dt oz oz 
dp" _ gd" 
az go” 
og! _ 36! 
oz) Cyan 


where ( ) denotes the time average operator. 
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(4.20) 


(4.21) 


(4.22) 


(4.23) 


(4.24a) 


(4.24b) 


Define the explicit terms at time / as 


ou l 
fees or (4.25) 
OV 
r, = Aok +fu+(v+ VA (4.26) 
OZ dy 
/ 
a, (4.27) 
oz 


Using (4.25)44.27) in (4.20)-(4.23), the following set of equations ts obtained: 





d,u 

H +r, - 0, (4.28) 
div do 
En, 00, - 0, (4.29) 
dt 
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Next. define the auxiliary variables 


is F(y,Z,t+ At), 
F° E F(y -a,z,t), 


F- =- F(y-2a,z,t-Ad, 
where 


a™! = At.(v°(y-a*,0+V)(y-a",d). 


(4.30) 


(4.31) 


(4.32a) 


(4.32b) 


(4.32c) 


(4.33) 


Note that the use of the time average of the vertical velocity w eliminates the 


need of considering vertical displacements in the calculation of the trajectories of the 


fluid particles. Therefore, the interpolation will be performed only tn the horizontal 
direction using one dimensional cubic spline functions. 
Define the following approximations tor the time derivatives and tame 


average respectively: 





dy? 2 D (4.34) 

dt 2At ` 

Fi (4.35) 
2 


Use of approximations (4.34) and (4.35) into (4.28)-(4.31) gives 


Ea mol (4.36) 
ŻAt 
TM ee (4.37) 


2At 2% ww) ' 








2 a (4.38) 














gy Oz 
wo w | 0, (4.38b) 
Oy az 
9-0 1 O (wi +w +r, a. (4.39) 
2At 2 & 


Group the terms at time (¢-At) and define new auxiliary variables p, ; 








b, see (4.40) 
j- 
P, = ~y +At dp ; (4.41) 
dy 
oe". 
= -0 +At w”. (4.42) 
P3 = 


28 


Similarly, group the terms at time /t + 4t)and define the auxiliary variables 


q ur, (4 43) 


j+ 
q - va 20 (4.44) 





$ 


00 ) 
- 0" +At w”. (4.45) 
q) 32 





Use of (4.40)-(4.45) allows to rewrite (4.36), (4.37) and (4.39) as follows: 


A ŻAR (4.46) 
qea-p,-ŻAtn, (4.47) 
Gq, > Pa 20. (4.48) 


The variables p, and r, can be calculated explicitly and with their values the 


q, variables can be obtained. 


To calculate the values of the variables u. v. w. 8 and q at the time level 
(¢ + At) it is necessary to solve the system of equations composed by (4.38b), (4.43)- 
(4.45) and the hydrostatic equation (4.24b). The procedure used by Robert et al.( 1985) 
consists of first solving an elliptic equation for q. 

In this model it ts not possible to state exact boundary conditions for g, due 
to the use of the rigid lid assumption. However, the w field has the exact boundary 
conditions (3.13). The following elliptic equation for w can be derived from the system 
of equations at time level (z + At): 


Sw A At? g 0" w* Atg oq, Mi q, (4.49) 


dz O & əy? do ay? aye 








Once w was calculated. 6 was determined using (4.45), and > and y were 
obtained using the procedure described in Chapter III. After adjusting the v field to 
satisfy < v > = 0, a new dp” field was calculated using the new values of v in (4.44) 
and the wfield was also adjusted to satisfy the continuity equation (4.38a). 

Finally, both in the explicit and the semi-implicit schemes the spatial 
derivatives are approximated by centered in space finite differences. Following Monk 
(1989), the derivatives are first evaluated at the grid points and, after that the values are 
interpolated for the position of the fluid particle. This procedure is computationally 
more efficient than the calculation of the derivatives at the fluid particle points, which 


would double the number of interpolations. 
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V. EXPERIMENTS AND RESULTS 


A. PARAMETERS USED IN THE EXPERIMENTS 


The following values for the constants were used in all experiments. 





L — 3600 km, 

H = 12 km. 

Ros ok. 

9,=300ºK, 

g=9.81ms 

OS 

DANOS 

z, = 9 km, 

aK fen’ 
Oz 


For the semi-Lagrangian semi-implicit (SLSI) case there were two control runs 
using Ay = 20 km and Az = 167 m. All the remaining experiments used normal 
resolution. with A y= 40 km and Az= 333m. The cross-sections of the initial fields are 
shown in Figs. 2-5. There is a warm front with thermally indirect circulation near y= 
900 km and a cold front with termally direct circulation near y= 2700 km. 

In order to assess the evolution of the frontogenesis process, the following 


parameter is defined at the lowest computational level: 


y - 148 


2 | 
Ya 
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Figure 2. Initial y component field (m/s). 


Figure 3. Initial vcomponent field (m/s). 
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Figure 4. Initial w component field (m/s). 





5. Initial 6 field (“K). 


Figure 


aS 


where A@ is the maximum variation ot 6 along the v direction. This parameter. used 
by Williams et al.(1991), gives a reasonable measure of the width of the frontal zone. 

During the first experiment wit he semi-implicit scheme. a high frequency 
oscillation was observed in the frontal width. That noise was not present in the 
experiments carried out with the semi-Lagrangian explicit (SLEX) scheme. In order to 
eliminate the noise, a time filtering technique similar to the one described by Asselin 


(1972) was used. The time filtering of a variable Q is defined as 


Q(t) — Q(D-vIQ(tAD-20(0-Q(t-AD), ĠA 


where Q(t+ 4 t)has been previously obtained by using the respective predictive equation 
with the unaveraged value of Q/t). The time filter defined above is expected to affect the 
higher frequencies only. 

Several experiments were performed with the objective of assessing how the use 


of time filtering would affect the solutions, as will be described next. 


B. DESCRIPTION OF THE EXPERIMENTS 
The experiments basically consisted of carrying out time integrations of the model 
using the SLEX and SLSI schemes. starting from the same initial conditions and using 
different combinations of time step A ¿and filtering factor y. The experiments were 
intended to give information about the following points: 
¢ How well the different schemes handle the formation of large gradients; 
e How the use of time filtering affects the final solutions; and 


e How the solutions are affected by the use of longer time steps. 


Table | shows the experiments performed with the SLEX scheme, asa function 
of time step Az and filtering factor y. The experiments carried out using the SLSI 


scheme are tdentified on Table 2. 


TABLE !. EXPERIMENTS WITH THE SLEX SCHEME. 


EM SLEX04 SLEXOS SLEX06 


TABLE 2. EXPERIMENTS WITH THE SLSI SCHEME. 


Joam Jus [same] = l- 

EN fase] — 

ssw fam [ans] 
Eq 
E 











p= - - 


SLSI16 | SLSII7 SLSI18 | SLSTI9 
SLSDO | SLSDI SLSD2 | SLSD6 | SLSL? 


Hel 


E 
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There were two high-resolution experiments. both using the SLSI scheme, with 
time steps of 180 secs. In the first of them , identified by SLHRO1, the time filtering 
factor y = 0.01 was used . and in the second one (SLHRO2) y was set equal to 0.05. The 
velocity components and temperature fields at time ¢ = 40 hours, for experiment 
SLHROI are shown tn Figs. 6-9 The occurrence of frontogenesis for the surface cold 
front and frontolysis for the warm front can be seen. These processes are due to the 
relation between the wind deformation and temperature disturbance fields, that gives 
convergence in the cold front region and divergence for the warm front. The wand v 
fields give a direct circulation around the cold frontal zone and an indirect circulation 
about the warm frontal zone. The u field develops cyclonic shear at the cold front. 

The evolution of the frontal width (d-va/ue) for SLHROL is presented in Fig. 10. 
An almost linear decrease in the frontal scale can be seen. This curve is different from 
that obtained by Williams et al.(1991), because in that case there was a momentum 
diffusion term that gave an asymptotic behavior in the evolution of the frontal width 
towards a balance condition. In the present model a continuous reduction of the frontal 
scale is expected until a limit, when the grid resolution is reached. Figure | | shows the 
evolution of the d-va/ue for experiment SLEXO1. It can be seen that the minimum value 
of the frontal width is reached at approximately ¢= 40 hours. Note that this curve ts 
more linear than that obtained in experiment SLHRO!. Also, the d-va/ue curves 
obtained for experiments SLEX02 and SLEX03 were virtually coincident with the one 
for SLEXO1. This result shows that the use of stronger filtering effect (larger y) does not 
affect the evolution of the frontogenesis process represented by the SLEX scheme. 

During the first experiment with the SLSI scheme (SLSIO1L) the d-va/ue curve. 
shown in Fig. 12, presented a high frequency oscillation. The time integration could be 


carried out until time ¢= 36 hours and after that, the iterative procedures no longer 
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Figure 6. u component field (m/s) for SLHROI at t — 40 hours. 





Figure 7. vcomponent field (m/s) for SLHROI att = 40 hours. 
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Figure 8. wcomponent field (m/s) for SLHROI at t = 40 hours. 
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Figure 9. 8 field (“K) for SLHRO1 at t = 40 hours. 
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Figure 10. Evolution of the frontal width (d-value) for SLHROI. 
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Figure 11. Evolution of the frontal width (d-value) for SLEXO1. 





39 


achieved convergence and the integration was terminated. It is important to note that 
the irregular behavior of the curve was not related to continuous amplification of the 
solution due to numerical instability and that the curve presents a general tendency to 
the frontogenesis solution obtained in the SLEX experiments. The high resolution 
experiment SLHRO1 also shows a low frequency oscillation about its linear decrease in 
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Figure 12. Evolution of the frontal width (d-value) for SLSIO1. 





In order to eliminate the noise the time filter described in the last section was 
introduced. Figure 13 shows the evolution of the frontal width obtained in experiment 
SLSIO2, with the filtering factor y set equal to 0.01. It can be seen from the figure that 
although the filtering effect was small, it was adequate to suppress the high frequency 


oscillations. 
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Figure 13. Evolution of the frontal width (d-value) for experiment SLSIO2. 


Figures 11 and 13 show that both for the the SLSI and SLEX experiments the 
minimum frontal width is reached approximately at time ¢= 40 hours and after that the 
d-value curves oscillate. Those oscillations can be related to the fact that the front 
reached a width corresponding to the model resolution and after that point the solutions 
were no longer physically consistent. 

Tables 3 and 4 show the minimum value of the frontal width for the SLEX 
experiments with time steps 180, 360 and 600 sec., and the corresponding SLSI 
experiments for which the time integrations could be carried out until /= 40 hours. 
Also presented is the time when the minimum d-va/ue ocurred. 

Comparison of the values presented in Tables 3 and 4 show that the minima d- 
values obtained in the SLSI experiments are smaller than those obtained for the 
corresponding SLEX experiments with same A fand y. Also, in the SLEX experiments 


the minimum frontal width increases as A fincrease whereas in the SLSI case there is no 
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significant change in the minimum d-va/ue due to changes in the time step. lt is 
interesting to note that for the SLEX case a shorter time integration is necessary for the 
minimum to occur for increasing time step. while the inverse tendency ts observed in the 
SLSI case that is. a longer time integration is necessary for the minimum frontal width 


to occur, for larger At 


TABLE 3. MINIMUM D-VALUE FOR SLEX EXPERIMENTS. 


A t (sec.) Experiment Minimum Time (hour) 
d-value (km) 


SLEXO1 
SLEXO2 
SLEX03 


SLEXOS 
SLEX06 
SLEXO7 


SLEXOS 


men 





f SA o a ee a ee 


Therefore, it can be concluded that the SLSI scheme gives a better representation 
of the scale collapse process than that obtained from the SLEX scheme, since narrower 


frontal widths can be resolved by the former than by the latter. 
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TABLE 4. MINIMUM D-VALUE FOR SLSI EXPERIMENTS. 
Minmum Time 


| At (sec.) Experiment 
d-value (km) | (hours) | 
180 CM SLSIO2 43.4 


a [as | 
ofa foe ur 
or fem fo [oe | 
o [sm fer je 
o fo for foe 


Experiment SLSIO3 used y= 0.03 and the d-va/uecurve obtained was practically 








coincident with the one obtained in SLSIO2, with y=0.01. 

In experiments SLSIO4 and SLSIO7 no filtering was used (y=0.0) and 
approximately at time ¢= 37 hours the numerical integration was interrupted because 
the solutions were no longer able to achieve convergence in the iterative procedures. 
However, the use of larger values of y (SLSIOS. SLSI06, SLSI24, SLSI08, SLSIO9. 
SLSI25) allowed the normal execution of the time integrations. Figure 14 shows the d- 
value curves for experiments SLSI08 and SLSIO9. [t can be seen that the oscillations 
were not completely filtered out with y=0.01 and a stronger filtering effect (y=0.03) was 
required. The experiments with longer time steps showed that as the time step increased 
a larger value for y was required for an appropriate filtering of the noise. Also, all 
experiments with the SLSI scheme showed that the largest amplitudes of the noise 


occurred during the first 24 hours of integration. Furthermore. those experiments 
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showed that the noise was not related to imbalance in the initial conditions, since the 
mass and wind fields were initially adjusted using the quasi-geostrophic circulation 
equation (4.8) and the oscillations were not present before f= 02 hours. The presence 
of the noise was investigated by examining the u, v and 6 fields obtained at the lowest 
computational level of experiment SLSIO7. At time ¢ =13 hours, that corresponds 
approximately to maximum amplitude of the noise, there was no evidence of small scale 
spatial oscillations in those fields. Although the origin of the noise could not be 
determined, the experiments showed that the use of the time filter was effective in 
eliminating the oscillations with relatively small values of y. It is expected that this 


filtering would not significantly affect the lower frequency solutions. 
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Figure 14. Evolution of the fronal width (d-value) for experiments SLSI08 and SLSIO9. 


Another point of interest in this study was the computational! effort employed in 
each of the experiments, since one of the objectives of the SLSI scheme is to allow stable 


solutions with long time steps. In order to evaluate the relative computational! 


efficiency, the CPU time spent for running each experiment for 42 hours was measured 
and normalized by the CPU time used in SLSIOL It was observed that the CPU time 
was not affected by changes in the value of y. so that the results are presented in [able 
5 asa function of numerical scheme and time step only. As expected, the computational 
effort decreases almost linearly for increasing time step. Also, for a given time step. the 
SLEX scheme ts more efficient than the SLSI one. This ts observed because the SI SI 
technique requires the solution of an elliptic equation every time step. [t should also be 
pointed out that in this study the SLSI formulation makes the spatial interpolations 
necessary only in the horizontal direction, whereas in the SLEX scheme the 
interpolations are performed using bicubic splines. Thus, a formulation of the SLSI 
scheme using two-dimensional interpolations could give a further increase in the 
computational effort. The advantage of the SLSI scheme appears clearly when 
integrations with time steps as long as 3600 sec. are performed with perfectly stable 
solutions, while the SLEX scheme did not allow time steps longer than 600 seconds. 

In order to assess the accuracy of the several solutions, a second high resolution 
(SLHRO2) control run was carried out with the SLSI scheme. The parameters used 
were At = 180 sec. and y = 0.05. This value for y was chosen to guarantee that the 
solutions would not be affected by noise. The solutions of the control run were linearly 
interpolated for the normal resolution grid point positions when necessary, because of 
the staggering of the variables. The differences between the values obtained for the 
control run and those obtained for the experiments with normal resolution were 
calculated at each grid point. The results were used to calculate the domain RMS 
differences for each variable. Two sets of comparisons were performed: In the firstone 
the velocity components and temperature fields for selected SLSI experiments at ¢= 36 


hours were compared with the ones obtained from SLHRO? at the same time. That 
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time was considered appropriate to give a well characterized atmospheric front and 
sufficiently away from the inconsistent solutions obtained after the model resolution 


was achieved. 


TABLE 5. NORMALIZED CPU TIME FOR 42-HOUR MODEL INTEGRATION. 


SLEX | SLSI 


18 


The experiments chosen for this first set of comparisons were those with y= 0.05 





for all time steps, except A 3600 sec., where y= 0.07 was used. The accuracy of the 
solutions was expressed in terms of RMS differences in the values the of u, v.w and 6, 
given in Table 6. 

The values obtained show that the u and 6 fields have small errors as compared 
to the range of the values observed in the domain (-35 to +22 °K for 6 and -50 to +15 
ms for u). On the other hand, the errors obtained in the v and w fields were relatively 


large compared with the magnitude of the reference values obtained for SLHRO2. The 
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large differences can be associated with the fact that the » and # components are 
closely related to the divergence in the model and their values can be significantly 
affected by the presence of gravity waves. Time series of y and + were examined and 
they confirmed the presence of tnternal gravity waves in the solutions. Thus, a single 
sample would not be sufficiently representative of the actual level of accuracy obtained 


in those fields. 
TABLE 6. RMS DIFFERENCES FOR TIME T= 36 HOURS. 
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SLSI23 and SLSTII 5 are shown in Figs. 15-18. It can be seen that the largest differences 
occur close to the frontal region for the u field, but .hey are almost uniformily 


distributed at the lowest levels in the 6 field for experiment SLSI15. 
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Figure 15. Differences in y component field (m/s) for experiment SLSI23 at t = 36 
hours. 





Figure 16. Differences in v component field (m/s) for experiment SLSIIS at t = 36 
hours. 
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Figure 17. Differences in 6 field (“K) for experiment SLSI23 at t = 36 hours. 





Figure 18. Differences in @ field ((K) for experiment SLSI15 at t = 36 hours. 
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The second set of comparisons used as reference the same solutions obtained for 
experiment SLHRO? at time ¢= 36 hours. The frontal width at that time in the high 
resolution experiment was 99 5 km. For those SLSI experiments where the time filter 
controlled the high frequency oscillations, the frontal width d= 100 km was uniquely 
related to a certain time. If the solutions had no error we could suppose that for a 
certain d-va/ue the velocity components and temperature fields should be the same as 
those obtained from the control run. Thus, the geophysical fields for the times 
corresponding to a d-va/ue approximately equal to 100 km were chosen for the second 
set of comparisons. Table 7 gives the RMS differences obtained and the corresponding 
frontal widths and times used for each comparison. The values obtained show that the 
frontogenesis process 1s slower for longer time steps. since the frontal width of 100 km 
is reached at a later time for larger Af. The RMS differences for u and 6 are again 
relatively small compared to the magnitude of the values obtained in the reference 
solution. The differences for v and w are still relatively large with respect to the 
reference values, as expected. Figures 19-22 show the grid point differences of the u 
and 6 fields for experiments SLSI23 and SLSIIS. 

The figures show that the largest differences occur close to the front. This 
characteristic of semi-Lagrangian schemes in which the largeerrors areconfined around 
the scale collapse region was observed by Kuo and Williams (1990). Such a property 
is desirable, since the regions away from the front will not undergo a large effect of 
errors generated close to the region where strong gradients are present. 

Figures 23-29 show the temperature distribution at the lowest computational level 


(z= 167 m) for the experiments listed in Table 5 (solid line) for comparison with the 
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values obtained from the control run (dashed line). [t can be seen that as the time step 
increases the differences increase but the largest errors remain close to the frontal 
region. 


TABLE 7. RMS DIFFERENCES FOR d ~ 100 km. 


o: RMS DIFFERENCES E 
Experiment, | u(m/s) | díkm) | Time 

| Atí(sec.), y a (hours) 
SLSIO3, 0.667 
180, 0.03 

I SLSIO6, 

| 360, 0.03 


| SLSIL?, 
| 1200, 0.05 
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Figure 19. Differences in y component field (m/s) for experiment SLSIO3 at time 
corresponding to d=100 km. 





Figure 20. Differences in u component field (m/s) for experiment SLSIIS at time 
corresponding to de100 km. 
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Figure 21. Differences in 6 field ("K) for experiment SLSI03 at time corresponding to 
d~100 km. 





Figure 22. Differences in 0 field (“K) for experiment SLSI15 at time corresponding to 
d~100 km. 
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For this model the Courant number a iscalculated using the following expression: 


© = Whar a al cad TDT) 5 (5.3) 
where C ma ÍS the phase speed of the fastest internal gravity wave. An important result 
from this study is that for the SLSI experiments, g ranged from 0.5 (A £= 180 sec.) to 
10.4 (A t= 3600 sec.) , as opposed to the CFL stabilty criterion that would require ø to 
be less than or equal to 1, and the values of the RMS differences obtained for u and 6 
in the experiments with large time steps remained relatively small. This result indicates 
that the gain in computational efficiency does not affect significantly the accuracy of the 


solutions obtained with the SLSI scheme. 
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Figure 23. Temperature disturbance at z = 167 m, for experiment SLSI03 (solid line) 
and control run (dashed line), at time corresponding to d=100 km. 
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Figure 24. Same as Fig. 23 except for experiment SLSIO6. 
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Figure 25. Same as Fig. 23 except for experiment SLSIO9. 
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---- Control Run 
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Figure 26. Same as fig. 23 except for experiment SLSI12. 
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Figure 27. Same as Fig. 23 except for experiment SLSIIS. 
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---- Control Run 
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Figure 28. Same as Fig. 23 except for experiment SLSI19. 
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Figure 29. Same as Fig. 23 except for experiment SLSI27. 
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VI. SUMMARY AND CONCLUSIONS 


In this study a numerical model based on the hydrostatic Boussinesq equations 
was used to simulate atmospheric frontogenesis driven by an irrotational non-divergent 
deformation wind field. The semi-Lagrangian technique was employed to integrate 
numerically the prognostic equations. The model assumes periodic boundary 
conditions in the cross-front direction and the rigid lid assumption is used as the upper 
boundary condition. The model also neglects along front variations. A basic sinusoidal 
deformation flow ts introduced as the forcing of the frontogenesis process. Experiments 
were performed using the semi-Lagrangian technique associated with two different time 
schemes: explicit and semi-implicit. In the semi-Lagrangian explicit case (SLEX) 
bicubic splines were used to interpolate the variables in space. In the semi-Lagrangian 
semi-implicit case (SLSI) the variables were interpolated in the y-direction only, using 
one-dimensional splines. A frontal width parameter (d-va/ue ) was defined at the lowest 
computational level to describe the evolution of the frontogenetical processes. The 
experiments with the SLEX technique were successful in replicating the formation of a 
realistic front in approximately 40 hours. Different time steps were used for the 
integration of the model and solutions which were both numerically and physically 
consistent were obtained with the SLEX scheme for values of Courant numbers as large 
as 1.7. 

The experiments with SLSI scheme presented high frequency oscillations that were 
eliminated by using a time filter. However, no small scale spatial oscillations were 


observed in the solutions. Different filtering effects were tested by changing the value 
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of the filtering parameter y. The solutions obtained with the SLEX scheme were not 
affected by the time filter, whereas the SLS[ solutions showed that stronger filtering was 
necessary for larger time steps. 

The SLSI scheme showed to be more successful in handling the scale collapse 
process than the SLEX scheme, since the fronts reproduced by the former had 
minimum widths smaller than those obtained by the latter. The SLSI scheme presented 
a tendency of slowing down the frontogenesis process for increasing time steps. 

Experiments were performed with the SLSI scheme with time steps as long as 3600 
seconds, corresponding to a Courant number greater than 10. The solutions obtained 
were perfectly stable, and the accuracy was not significantly degraded, even for longer 
time steps. 

The SLSI solutions also had the characteristic of constraining the largest errors 
close to the scale collapse region. Such a characteristic is desirable since the errors will 
have a smaller impact on the solutions in regions away from the front. 

The results obtained in this study suggest that the SLSI technique is appropriate 
for mesoscale regional! models since it is computationally efficient and produces accurate 
results. A different formulation of the scheme where two-dimensional interpolation of 
variables were allowed should be studied to verify if there would be a significant positive 
impact on the accuracy of the solutions. Also, effects of surface topography and 
physical processes like advection of the frontal system, friction, vertical shear of the 


basic wind field and moisture should be investigated in future studies. 
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